SEASONAL TEMPERATURES

Figure 1a: Temperatures are seasonal

temp3=read.csv("foweytemps5yrs.csv")
temp3$date=as.Date(with(temp3, paste('2019',MM, DD,sep="-")), "%Y-%m-%d")
 temp3$WTMP=as.numeric(paste(temp3$WTMP))
 temp=read.csv("temp.csv")
temp$Date=as.Date(temp$Date,format = "%d/%m/%Y")
 
 ggplot()+stat_summary(data=temp3,aes(x=date,y=WTMP,fill='darkgreen'),fun.data=mean_se,geom='ribbon',alpha=0.6)+theme_minimal(base_size = 15)%+replace%
    theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank())+
   labs(x='Month (2019-20)',y='Temperature/°C')+stat_summary(data=temp,aes(x=Date,y=Temp,fill='chartreuse2'),fun.data=mean_se,geom='ribbon',alpha=0.6)+scale_fill_identity(guide='legend',breaks=c('darkgreen','chartreuse2'),labels=c('Fowey Rock 2015-19 average','Emerald Reef 2019'))+guides(fill=guide_legend(title=''))+theme(legend.position = c(0.6,0.3))+scale_x_date(date_labels = '%b',limits=as.Date(c('2019-01-01','2019-12-31')))

#ggsave('temp3.pdf',device='pdf', width=7,height=5) #add labels to figure to show two collection timepoints

INITIAL SYMBIONT COMMUNITIES

Figure 1b: M. cavernosa symbionts per host cell increased from April to October.

Add arrows to show that points above line show higher in april, below line shows higher in April.

#calculate seasonal per colony change in symbiont density (ie the average symbiont density between the two cores taken from each colony per batch)
batches$pre_sh[which(batches$pre_sh==0)]=NA
Apr_Oct_sh= batches %>%
  group_by(Colony,Batch) %>%
  summarise_at(vars(pre_sh), funs(mean(., na.rm=TRUE)))

Apr_sh= filter(Apr_Oct_sh, Batch=='April')
Apr_sh$Apr_pre_sh=paste(Apr_sh$pre_sh)
Oct_sh= filter(Apr_Oct_sh, Batch=='October')
Oct_sh$Oct_pre_sh=paste(Oct_sh$pre_sh)
Apr_Oct_sh_change=join(Apr_sh,Oct_sh,by='Colony')
Apr_Oct_sh_change=Apr_Oct_sh_change[,c(1,4,7)]
Apr_Oct_sh_change$Apr_pre_sh=as.numeric(Apr_Oct_sh_change$Apr_pre_sh)
Apr_Oct_sh_change$Oct_pre_sh=as.numeric(Apr_Oct_sh_change$Oct_pre_sh)
Apr_Oct_sh_change$pre_sh_change= Apr_Oct_sh_change$Oct_pre_sh- Apr_Oct_sh_change$Apr_pre_sh
Apr_Oct_sh_change$rel_sh_change= ((Apr_Oct_sh_change$Oct_pre_sh- Apr_Oct_sh_change$Apr_pre_sh)/Apr_Oct_sh_change$Apr_pre_sh)
Apr_Oct_sh_change$octtoapr<- Apr_Oct_sh_change$Oct_pre_sh/Apr_Oct_sh_change$Apr_pre_sh
Apr_Oct_sh_change$Species=factor(Apr_Oct_sh_change$Colony)
Apr_Oct_sh_change$Species= mapvalues(Apr_Oct_sh_change$Species,
                          from=c('100','2','27','28','39','67','68','71','72','87'),to=(rep('M.cavernosa',times=10)))
 Apr_Oct_sh_change$Species= mapvalues(Apr_Oct_sh_change$Species,
                          from=c('13','21','34','36','4','62','64','65','66','81'),to=(rep('O.faveolata',times=10)))                         
  Apr_Oct_sh_change$Species= mapvalues(Apr_Oct_sh_change$Species,
                          from=c('16','18','19','22','23','26','3','35','41','48'),to=(rep('S.siderea',times=10))) 
# this dataframe shows the average change in s:h per colony


##test statistical significance in seasonal change in S:H
#hist(log10(batches$pre_sh)) #log10 transformed data then used linear mixed effects model on s:h data to test effect of batch within each species, with colony as a random factor
batches$transf_presh= log10(batches$pre_sh) # transform the response variable

mcavpreshmod=lmer(log10(pre_sh)~Batch+(1|Colony),data=filter(batches,batches$Species=='M.cavernosa'))
plot_resqq(mcavpreshmod) # check residuals are normally distributed

mcavpreshmod2<- as_lmerModLmerTest(mcavpreshmod)
summary(mcavpreshmod2) #p for batch, p=0.001126** 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(pre_sh) ~ Batch + (1 | Colony)
##    Data: filter(batches, batches$Species == "M.cavernosa")
## 
## REML criterion at convergence: 56
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.17638 -0.65717 -0.08901  0.55392  2.38020 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Colony   (Intercept) 3.240e-08 0.00018 
##  Residual             2.357e-01 0.48546 
## Number of obs: 38, groups:  Colony, 10
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   -2.3177     0.1144 35.9720  -20.25  < 2e-16 ***
## BatchOctober   0.5583     0.1577 35.9823    3.54  0.00113 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BatchOctobr -0.725
# now let's get the batch parameter estimates and CIs:
mcavemm.sh<- emmeans(mcavpreshmod, specs=revpairwise~Batch) 
summary(mcavemm.sh)
## $emmeans
##  Batch   emmean    SE   df lower.CL upper.CL
##  April    -2.32 0.115 25.7    -2.55    -2.08
##  October  -1.76 0.109 25.7    -1.98    -1.54
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast        estimate    SE   df t.ratio p.value
##  October - April    0.558 0.158 28.7 3.525   0.0014 
## 
## Note: contrasts are still on the log10 scale 
## Degrees-of-freedom method: kenward-roger
mcavemmsh_contrasts<- mcavemm.sh$contrasts %>%
     confint()%>%
     as.data.frame() # NB these results are given on a log10 scale, giving the odds ratio of October to April S:H

mcavemm.sh2<- emmeans(mcavpreshmod, specs=revpairwise~Batch, type='response') 
summary(mcavemm.sh2) #alternatively, by backtransforming the predicted means for each batch, we can then calculate the percentage change, which may be more intuitive to interpret.
## $emmeans
##  Batch   response      SE   df lower.CL upper.CL
##  April    0.00481 0.00128 25.7  0.00279  0.00831
##  October  0.01741 0.00435 25.7  0.01041  0.02910
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log10 scale 
## 
## $contrasts
##  contrast        estimate    SE   df t.ratio p.value
##  October - April    0.558 0.158 28.7 3.525   0.0014 
## 
## Note: contrasts are still on the log10 scale 
## Degrees-of-freedom method: kenward-roger
mcavemmsh_pred<-as.data.frame(mcavemm.sh2$emmeans)

mcavemmsh_pred$relchange<- (mcavemmsh_pred[2,2]-mcavemmsh_pred[1,2])/mcavemmsh_pred[1,2]
mcavemmsh_pred$relupper<- (mcavemmsh_pred[2,6]-mcavemmsh_pred[1,5])/mcavemmsh_pred[1,2]
mcavemmsh_pred$rellower<- (mcavemmsh_pred[2,5]-mcavemmsh_pred[1,6])/mcavemmsh_pred[1,2]



ofavpreshmod=lmer(log10(pre_sh)~Batch+(1|Colony),data=filter(batches,batches$Species=='O.faveolata'))
plot_resqq(ofavpreshmod) # check residuals are normally distributed

ofavpreshmod2<- as_lmerModLmerTest(ofavpreshmod)
summary(ofavpreshmod2) #p for batch, p=0.938
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(pre_sh) ~ Batch + (1 | Colony)
##    Data: filter(batches, batches$Species == "O.faveolata")
## 
## REML criterion at convergence: 29.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3754 -0.5262  0.1305  0.6129  1.6623 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.12853  0.3585  
##  Residual             0.07002  0.2646  
## Number of obs: 35, groups:  Colony, 10
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -1.491446   0.128795 11.205854 -11.580  1.4e-07 ***
## BatchOctober  0.007186   0.091378 24.399854   0.079    0.938    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BatchOctobr -0.320
# now let's get the batch parameter estimates and CIs:
ofavemm.sh<- emmeans(ofavpreshmod, specs=revpairwise~Batch) 
summary(ofavemm.sh)
## $emmeans
##  Batch   emmean    SE   df lower.CL upper.CL
##  April    -1.49 0.129 11.1    -1.77    -1.21
##  October  -1.48 0.132 12.1    -1.77    -1.20
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast        estimate     SE   df t.ratio p.value
##  October - April  0.00719 0.0915 24.3 0.078   0.9381 
## 
## Note: contrasts are still on the log10 scale 
## Degrees-of-freedom method: kenward-roger
ofavemmsh_contrasts<- ofavemm.sh$contrasts %>%
     confint()%>%
     as.data.frame() # NB these results are given on a log10 scale

ofavemm.sh2<- emmeans(ofavpreshmod, specs=revpairwise~Batch, type='response') 
summary(ofavemm.sh2) #alternatively, by backtransforming the predicted means for each batch, we can then calculate the percentage change, which may be more intuitive to interpret.
## $emmeans
##  Batch   response      SE   df lower.CL upper.CL
##  April     0.0323 0.00957 11.1   0.0168   0.0619
##  October   0.0328 0.00997 12.1   0.0169   0.0635
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log10 scale 
## 
## $contrasts
##  contrast        estimate     SE   df t.ratio p.value
##  October - April  0.00719 0.0915 24.3 0.078   0.9381 
## 
## Note: contrasts are still on the log10 scale 
## Degrees-of-freedom method: kenward-roger
ofavemmsh_pred<-as.data.frame(ofavemm.sh2$emmeans)

ofavemmsh_pred$relchange<- (ofavemmsh_pred[2,2]-ofavemmsh_pred[1,2])/ofavemmsh_pred[1,2]
ofavemmsh_pred$relupper<- (ofavemmsh_pred[2,6]-ofavemmsh_pred[1,5])/ofavemmsh_pred[1,2]
ofavemmsh_pred$rellower<- (ofavemmsh_pred[2,5]-ofavemmsh_pred[1,6])/ofavemmsh_pred[1,2]


ssidpreshmod=lmer(log10(pre_sh)~Batch+(1|Colony),data=filter(batches,batches$Species=='S.siderea'))
plot_resqq(mcavpreshmod) # check residuals are normally distributed

ssidpreshmod2<- as_lmerModLmerTest(ssidpreshmod)
summary(ssidpreshmod2) #p for batch, p=0.04833*
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(pre_sh) ~ Batch + (1 | Colony)
##    Data: filter(batches, batches$Species == "S.siderea")
## 
## REML criterion at convergence: 55
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.53246 -0.34681  0.00358  0.49889  1.82110 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.2705   0.5201  
##  Residual             0.1524   0.3903  
## Number of obs: 36, groups:  Colony, 9
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   -3.0287     0.1963 10.0494 -15.431 2.51e-08 ***
## BatchOctober  -0.2696     0.1301 26.0000  -2.072   0.0483 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BatchOctobr -0.331
# now let's get the batch parameter estimates and CIs:
ssidemm.sh<- emmeans(ssidpreshmod, specs=revpairwise~Batch) 
summary(ssidemm.sh)
## $emmeans
##  Batch   emmean    SE   df lower.CL upper.CL
##  April    -3.03 0.196 10.1    -3.47    -2.59
##  October  -3.30 0.196 10.1    -3.74    -2.86
## 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log10 (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast        estimate   SE df t.ratio p.value
##  October - April    -0.27 0.13 26 -2.072  0.0483 
## 
## Note: contrasts are still on the log10 scale 
## Degrees-of-freedom method: kenward-roger
ssidemmsh_contrasts<- ssidemm.sh$contrasts %>%
     confint()%>%
     as.data.frame() # NB these results are given on a log10 scale

ssidemm.sh2<- emmeans(ssidpreshmod, specs=revpairwise~Batch, type='response') 
summary(ssidemm.sh2) #alternatively, by backtransforming the predicted means for each batch, we can then calculate the percentage change, which may be more intuitive to interpret. However, this is misleading since it is the differnece of the averages, not the average difference. 
## $emmeans
##  Batch   response       SE   df lower.CL upper.CL
##  April   0.000936 0.000423 10.1 0.000342  0.00256
##  October 0.000503 0.000227 10.1 0.000184  0.00138
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log10 scale 
## 
## $contrasts
##  contrast        estimate   SE df t.ratio p.value
##  October - April    -0.27 0.13 26 -2.072  0.0483 
## 
## Note: contrasts are still on the log10 scale 
## Degrees-of-freedom method: kenward-roger
ssidemmsh_pred<-as.data.frame(ssidemm.sh2$emmeans)

ssidemmsh_pred$relchange<- (ssidemmsh_pred[2,2]-ssidemmsh_pred[1,2])/ssidemmsh_pred[1,2]
ssidemmsh_pred$relupper<- (ssidemmsh_pred[2,6]-ssidemmsh_pred[1,5])/ssidemmsh_pred[1,2]
ssidemmsh_pred$rellower<- (ssidemmsh_pred[2,5]-ssidemmsh_pred[1,6])/ssidemmsh_pred[1,2]

#collate the dataframe for the log10 transformed contrast ratios
sh_emmcontrasts<-rbind(mcavemmsh_contrasts,ofavemmsh_contrasts,ssidemmsh_contrasts)
sh_emmcontrasts$Species<-c('M.cavernosa','O.faveolata','S.siderea')
sh_emmcontrasts$Species<-as.factor(sh_emmcontrasts$Species)
sh_emmcontrasts$change_untransformed= 10^(sh_emmcontrasts$estimate)
sh_emmcontrasts$lowerCI_change_untransformed= 10^(sh_emmcontrasts$lower.CL)
sh_emmcontrasts$upperCI_change_untransformed= 10^(sh_emmcontrasts$upper.CL) #these ratios of oct:apr are now backtransformed and ready to be plotted 

#for the alternative plot, collate the calculated relative changes
sh_emmpredchange<-rbind(mcavemmsh_pred[1,],ofavemmsh_pred[1,],ssidemmsh_pred[1,])
sh_emmpredchange$Species<-c('M.cavernosa','O.faveolata','S.siderea')
sh_emmpredchange$Species<- as.factor(sh_emmpredchange$Species)


#plot set-up
ofav=expression(paste(italic("O. faveolata")))
ssid=expression(paste(italic("S. siderea")))
mcav=expression(paste(italic("M. cavernosa")))


ggplot()+
  geom_point(data=sh_emmcontrasts, aes(x=Species,y=change_untransformed, colour=Species), size=2)+ #large points are showing estimated october:april 
  geom_point(data=Apr_Oct_sh_change, aes(x=Species,y=octtoapr, colour=Species), size=0.5,  position=position_jitter(width=0.1))+ #small points are showing october s:h / april s:h
  geom_errorbar(data=sh_emmcontrasts, aes(x=Species, ymin=(lowerCI_change_untransformed), ymax=(upperCI_change_untransformed), colour=Species), size=0.5, width=0.2)+
  theme_minimal(base_size = 15)%+replace%
    theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank())+
  labs(y='October : April symbionts per host cell', x='')+
   scale_x_discrete(labels=c(mcav,ofav,ssid), expand=expansion(add=c(0.4,1.2)))+
  scale_color_manual(values=c('deeppink2','darkorange1', 'darkturquoise'))+
  guides(colour=F)+
  geom_hline(yintercept=0,linetype='dashed')

#ggsave('initialSH_seasonal_emmeans.pdf',device='pdf',width=7,height=5) # add 'higher in Oct'/'higher in Apr' labels to figure


#perhaps plotting on a log10 scale will correct the +ve -ve skew?
ggplot()+
  geom_point(data=sh_emmcontrasts, aes(x=Species,y=(estimate), colour=Species), size=2)+ #large points are showing estimated october:april 
  geom_point(data=Apr_Oct_sh_change, aes(x=Species,y=log10(octtoapr), colour=Species), size=0.5,  position=position_jitter(width=0.1))+ #small points are showing october s:h / april s:h
  geom_errorbar(data=sh_emmcontrasts, aes(x=Species, ymin=(lower.CL), ymax=(upper.CL), colour=Species), size=0.5, width=0.2)+
  theme_minimal(base_size = 15)%+replace%
    theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank())+
  labs(y='log(October : April symbionts per host cell)', x='')+
   scale_x_discrete(labels=c(mcav,ofav,ssid), expand=expansion(add=c(0.4,1.2)))+
  scale_color_manual(values=c('deeppink2','darkorange1', 'darkturquoise'))+
  guides(colour=F)+
  geom_hline(yintercept=0,linetype='dashed')

  ggplot()+
  geom_point(data=sh_emmpredchange, aes(x=Species,y=relchange, colour=Species), size=2)+ 
  geom_point(data=Apr_Oct_sh_change, aes(x=Species,y=rel_sh_change, colour=Species), size=0.5,  position=position_jitter(width=0.1))+ #points are showing proportional change from apr to oct
  geom_errorbar(data=sh_emmpredchange, aes(x=Species, ymin=rellower, ymax=relupper, colour=Species), size=0.5, width=0.2)+
  theme_minimal(base_size = 15)%+replace%
    theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank())+
  labs(y='relative change in symbionts per host cell', x='')+
   scale_x_discrete(labels=c(mcav,ofav,ssid), expand=expansion(add=c(0.4,1.2)))+
  scale_color_manual(values=c('deeppink2','darkorange1', 'darkturquoise'))+
  guides(colour=F)+
  geom_hline(yintercept=0,linetype='dashed')

#UPDATES FOR ROSS: 
#potential problem with this plot is that by plotting the ratio, a ratio of 0.5 (50% decrease) may be equivalent to a ratio of 2 (100% increase)...? Is this plot disproportionately inflating positive values? I have also now removed the scale limits so that every data point is shown.

Figure 1c: Proportion Durusdinium did not change between April and October for any of the 3 coral species.

batch_comparevi=
ggplot()+
   theme_minimal(base_size = 15)%+replace%
    theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank())+ 
  labs(y='Initial proportion *Durusdinium*')+
  scale_fill_manual(values=c('#3C5488B2','#00A087B2'),labels=c('April','October'))+
  geom_violin(data=batches,aes(x=Species,y=pre_propD,colour=Batch),scale = 'width')+
  geom_dotplot(data=batches,bins=30,binaxis='y',dotsize=0.7,stackratio=0.5,stackdir='center',stackgroups=F, position='dodge',aes(x= Species,y=pre_propD,fill=Batch))+
  scale_colour_manual(values=c('#3C5488B2','#00A087B2'),labels=c('April','October'))+
  scale_x_discrete(labels=c(mcav,ofav,ssid))+
   theme(axis.title.y.left  = element_markdown(), legend.position = c(0.2,0.5))+
  labs(x='')+
  theme(legend.title = element_text(size=0))
batch_comparevi

#ggsave('batchcompare_propd_vi_dot.pdf',device='pdf',width=7,height=5)

hist(batches$pre_propD)

mcavpredmod=glmer(pre_propD~Batch+(1|Colony),data=filter(batches,batches$Species=='M.cavernosa'), family = 'binomial')
summary(mcavpredmod) # p=1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: pre_propD ~ Batch + (1 | Colony)
##    Data: filter(batches, batches$Species == "M.cavernosa")
## 
##      AIC      BIC   logLik deviance df.resid 
##      6.0     11.1      0.0      0.0       37 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
##       0       0       0       0 2904892 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Colony (Intercept) 0        0       
## Number of obs: 40, groups:  Colony, 10
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -3.806e+01  1.501e+07       0        1
## BatchOctober -3.033e+02  2.122e+07       0        1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BatchOctobr -0.707
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
ofavpredmod=glmer(pre_propD~Batch+(1|Colony),data=filter(batches,batches$Species=='O.faveolata'), family = 'binomial')
summary(ofavpredmod) # p=0.2
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: pre_propD ~ Batch + (1 | Colony)
##    Data: filter(batches, batches$Species == "O.faveolata")
## 
##      AIC      BIC   logLik deviance df.resid 
##     31.9     36.5    -12.9     25.9       32 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.76647 -0.04036 -0.00779  0.25178  1.00309 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Colony (Intercept) 80.35    8.964   
## Number of obs: 35, groups:  Colony, 10
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -6.153      3.793  -1.622    0.105
## BatchOctober   -3.286      2.769  -1.187    0.235
## 
## Correlation of Fixed Effects:
##             (Intr)
## BatchOctobr 0.343
ssidpredmod=glmer(pre_propD~Batch+(1|Colony),data=filter(batches,batches$Species=='S.siderea'), family = 'binomial')
summary(ssidpredmod)# p0.09
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: pre_propD ~ Batch + (1 | Colony)
##    Data: filter(batches, batches$Species == "S.siderea")
## 
##      AIC      BIC   logLik deviance df.resid 
##     36.2     41.0    -15.1     30.2       33 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9098 -0.3598 -0.0362  0.2155  0.9639 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Colony (Intercept) 24.15    4.914   
## Number of obs: 36, groups:  Colony, 9
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    -1.039      2.132  -0.487   0.6259  
## BatchOctober    4.593      2.713   1.693   0.0905 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BatchOctobr -0.425
#Make it clear that this plot is showing raw data, not predictive model. glmer is now 'family quasibinomial'. 

BLEACHING SENSITIVITY & RATE

Figure 1d: M. cavernosa showed higher sensitivity to bleaching in October compared to April.

blch_sensitivity=blch_sensitivity[-c(56,57,75),] #remove 're'_drop_sh' NAs
mcavsensitivity=subset(blch_sensitivity,Species=='M.cavernosa')

mcavblchresmod=glmer(post_sh~Batch+(1|Colony),data=mcavsensitivity)
plot_resqq(mcavblchresmod)

summary(mcavblchresmod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: post_sh ~ Batch + (1 | Colony)
##    Data: mcavsensitivity
## 
## REML criterion at convergence: -171.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2063 -0.2934 -0.0358  0.0789  3.5715 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Colony   (Intercept) 3.308e-06 0.001819
##  Residual             4.659e-05 0.006826
## Number of obs: 27, groups:  Colony, 10
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.009471   0.002080   4.552
## Batch2      -0.009326   0.002664  -3.500
## 
## Correlation of Fixed Effects:
##        (Intr)
## Batch2 -0.718
mcavblchresmod2<- as_lmerModLmerTest(mcavblchresmod)
summary(mcavblchresmod2)#the effect of batch on post heat stress y2 is significant (0.002206).
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: post_sh ~ Batch + (1 | Colony)
##    Data: mcavsensitivity
## 
## REML criterion at convergence: -171.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2063 -0.2934 -0.0358  0.0789  3.5715 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Colony   (Intercept) 3.308e-06 0.001819
##  Residual             4.659e-05 0.006826
## Number of obs: 27, groups:  Colony, 10
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  0.009471   0.002080 19.873385   4.552 0.000196 ***
## Batch2      -0.009326   0.002664 20.380580  -3.500 0.002206 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## Batch2 -0.718
ofavsensitivity<-subset(blch_sensitivity,Species=='O.faveolata')
ofavblchresmod=glmer(post_sh~Batch+InitialDom+(1|Colony),data=ofavsensitivity)
plot_resqq(ofavblchresmod)

summary(ofavblchresmod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: post_sh ~ Batch + InitialDom + (1 | Colony)
##    Data: ofavsensitivity
## 
## REML criterion at convergence: -204.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.26588 -0.42748 -0.26007  0.03075  1.97818 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Colony   (Intercept) 9.239e-07 0.0009612
##  Residual             1.879e-06 0.0013706
## Number of obs: 24, groups:  Colony, 10
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)    1.286e-03  8.000e-04   1.608
## Batch2         2.554e-04  6.088e-04   0.420
## InitialDomnond 6.016e-05  8.994e-04   0.067
## 
## Correlation of Fixed Effects:
##             (Intr) Batch2
## Batch2      -0.380       
## InitilDmnnd -0.791  0.081
ofavblchresmod2<- as_lmerModLmerTest(ofavblchresmod)
summary(ofavblchresmod2) #when separated by symbiont genus, the effect of batch on post heat stress y2 is insignificant (0.680)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: post_sh ~ Batch + InitialDom + (1 | Colony)
##    Data: ofavsensitivity
## 
## REML criterion at convergence: -204.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.26588 -0.42748 -0.26007  0.03075  1.97818 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Colony   (Intercept) 9.239e-07 0.0009612
##  Residual             1.879e-06 0.0013706
## Number of obs: 24, groups:  Colony, 10
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)    1.286e-03  8.000e-04 8.141e+00   1.608    0.146
## Batch2         2.554e-04  6.088e-04 1.615e+01   0.420    0.680
## InitialDomnond 6.016e-05  8.994e-04 6.479e+00   0.067    0.949
## 
## Correlation of Fixed Effects:
##             (Intr) Batch2
## Batch2      -0.380       
## InitilDmnnd -0.791  0.081
ssidsensitivity<-subset(blch_sensitivity,Species=='S.siderea')
ssidblchresmod=glmer(post_sh~Batch+InitialDom+(1|Colony),data=ssidsensitivity)
plot_resqq(ssidblchresmod)

summary(ssidblchresmod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: post_sh ~ Batch + InitialDom + (1 | Colony)
##    Data: ssidsensitivity
## 
## REML criterion at convergence: -255.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6829 -0.4914 -0.0217  0.1180  2.4877 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Colony   (Intercept) 2.781e-08 0.0001667
##  Residual             4.270e-08 0.0002066
## Number of obs: 22, groups:  Colony, 9
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)     9.935e-05  1.135e-04   0.875
## Batch2          1.182e-04  1.002e-04   1.179
## InitialDomnond -7.638e-05  1.421e-04  -0.538
## 
## Correlation of Fixed Effects:
##             (Intr) Batch2
## Batch2      -0.579       
## InitilDmnnd -0.688  0.371
ssidblchresmod2<- as_lmerModLmerTest(ssidblchresmod)
summary(ssidblchresmod2) #when separated by symbiont genus, the effect of batch on post heat stress y2 is insignificant (0.258)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: post_sh ~ Batch + InitialDom + (1 | Colony)
##    Data: ssidsensitivity
## 
## REML criterion at convergence: -255.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6829 -0.4914 -0.0217  0.1180  2.4877 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Colony   (Intercept) 2.781e-08 0.0001667
##  Residual             4.270e-08 0.0002066
## Number of obs: 22, groups:  Colony, 9
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)     9.935e-05  1.135e-04  1.325e+01   0.875    0.397
## Batch2          1.182e-04  1.002e-04  1.414e+01   1.179    0.258
## InitialDomnond -7.638e-05  1.421e-04  1.330e+01  -0.538    0.600
## 
## Correlation of Fixed Effects:
##             (Intr) Batch2
## Batch2      -0.579       
## InitilDmnnd -0.688  0.371
mcavsensitivity$predicted<- predict(mcavblchresmod)
mcavsensitivity$residuals<- residuals(mcavblchresmod)
mcavsensitivityse<-mcavsensitivity %>%
  group_by(Batch) %>%
  summarise(se=std.error(post_sh))
as.data.frame(mcavsensitivityse)
##   Batch           se
## 1     1 3.069162e-03
## 2     2 5.373349e-05
mcavsensitivitysumm<-mcavsensitivity %>%
  group_by(Batch) %>%
  summarise(mean=mean(post_sh))
as.data.frame(mcavsensitivitysumm)
##   Batch         mean
## 1     1 0.0095320145
## 2     2 0.0001428238
mcavsensitivitysumm$se<-mcavsensitivityse$se
mcavsensitivitysey2<-mcavsensitivity %>%
  group_by(Batch) %>%
  summarise(sey2=std.error(post_y2))
as.data.frame(mcavsensitivitysey2)
##   Batch       sey2
## 1     1 0.01253447
## 2     2 0.01699856
mcavsensitivitymeany2<-mcavsensitivity %>%
  group_by(Batch) %>%
  summarise(meany2=mean(post_y2))
as.data.frame(mcavsensitivitymeany2)
##   Batch    meany2
## 1     1 0.1670833
## 2     2 0.2678667
mcavsensitivitysumm$sey2<-mcavsensitivitysey2$sey2
mcavsensitivitysumm$meany2<-mcavsensitivitymeany2$meany2


#UPDATES FOR ROSS: 
#what if we try plotting end fv/fm and s:h rather than the relative change? And use prediction ellipses or mean and error bars rather than try to fit a linear model. 

ggplot()+
  geom_point(data=mcavsensitivity,aes(x=post_y2,y=post_sh,colour=Batch))+
  geom_errorbar(data=mcavsensitivitysumm, aes(x=meany2,ymin=mean-se,ymax=mean+se, colour=Batch), width=0.01, show.legend=F)+
   geom_errorbar(data=mcavsensitivitysumm, aes(y=mean,xmin=meany2-sey2,xmax=meany2+sey2, colour=Batch), width=0.0025, show.legend=F)+
  geom_point(data=mcavsensitivitysumm, aes(x=meany2,y=mean,colour=Batch),size=3, show.legend=F)+
  theme_minimal(base_size = 13)%+replace%
    theme(panel.border = element_rect(colour='grey20',size=0.5,fill=NA))+
   scale_colour_manual(values=c('#3C5488B2','#00A087B2'),labels=c('April','October'))+
 theme(panel.spacing = unit(1.2, "lines"),legend.position = 'right')+
  theme(legend.title=element_text(size=0), axis.title.y = element_markdown())+
theme(legend.position = c(0.8,0.5))+
  labs(x='Fv/Fm after heat stress',y='symbionts per *M. cavernosa* cell after heat stress')

  #ggsave('bleaching_sensitivity_batches.pdf',device='pdf',width=7,height=5)
allbleach=rbind(ofav_DHW,ssid_DHW,mcav_DHW)
allbleach$Species=factor(allbleach$Species,levels=c('M.cavernosa','O.faveolata','S.siderea'))
allrecov=rbind(mcav_recov, ofav_recov, ssid_recov)
allrecov$Species=factor(allrecov$Species,levels=c('M.cavernosa','O.faveolata','S.siderea'))
allrecov$InitialDom=revalue(allrecov$InitialDom,c('D'='d','B'='nond','C'='nond'))
allrecov$InitialDom=factor(allrecov$InitialDom)
allbleach$InitialDom=revalue(allbleach$InitialDom,c('D'='d','B'='nond','C'='nond'))
allbleach$InitialDom=factor(allbleach$InitialDom)
allbleach$Batch=factor(allbleach$Batch)
allrecov$Batch=factor(allrecov$Batch)
allbleach$Colony=factor(allbleach$Colony)
allrecov$Colony=factor(allrecov$Colony)


allbleach<-allbleach[-c(167),] 
allrecov<-allrecov[-c(86,239),]#model highlighted this one datapoint as an outlier in both dataframes
  

bleachingmod2<-glmer(Shprop~DHW*Species*InitialDom+(1|Colony), data=allbleach)
plot_resqq(bleachingmod2) #perform statistical tests on mixed effects model

bleachingmod3<-as_lmerModLmerTest(bleachingmod2)
  summary(bleachingmod3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ DHW * Species * InitialDom + (1 | Colony)
##    Data: allbleach
## 
## REML criterion at convergence: 62.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1750 -0.4013 -0.0358  0.3297  5.1995 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.006468 0.08042 
##  Residual             0.060175 0.24531 
## Number of obs: 147, groups:  Colony, 29
## 
## Fixed effects:
##                                         Estimate Std. Error         df t value
## (Intercept)                             1.054054   0.125491  77.511446   8.399
## DHW                                    -0.059963   0.017904 120.365702  -3.349
## SpeciesO.faveolata                     -0.057977   0.159516  67.688528  -0.363
## SpeciesS.siderea                       -0.099000   0.102079  67.803938  -0.970
## InitialDomnond                         -0.023200   0.113848  86.620959  -0.204
## DHW:SpeciesO.faveolata                 -0.009859   0.020228 119.185291  -0.487
## DHW:SpeciesS.siderea                    0.002114   0.016147 121.484396   0.131
## DHW:InitialDomnond                     -0.044485   0.015250 121.587128  -2.917
## SpeciesO.faveolata:InitialDomnond      -0.033899   0.164764  68.629139  -0.206
## DHW:SpeciesO.faveolata:InitialDomnond   0.023099   0.019905 120.391812   1.161
##                                       Pr(>|t|)    
## (Intercept)                           1.65e-12 ***
## DHW                                    0.00108 ** 
## SpeciesO.faveolata                     0.71740    
## SpeciesS.siderea                       0.33557    
## InitialDomnond                         0.83901    
## DHW:SpeciesO.faveolata                 0.62690    
## DHW:SpeciesS.siderea                   0.89605    
## DHW:InitialDomnond                     0.00421 ** 
## SpeciesO.faveolata:InitialDomnond      0.83760    
## DHW:SpeciesO.faveolata:InitialDomnond  0.24814    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DHW    SpcsO. SpcsS. IntlDm DHW:SpO. DHW:SS DHW:ID SO.:ID
## DHW         -0.592                                                          
## SpecsO.fvlt -0.787  0.466                                                   
## SpecisS.sdr -0.798  0.536  0.627                                            
## InitilDmnnd -0.907  0.507  0.714  0.639                                     
## DHW:SpcsO.f  0.524 -0.885 -0.590 -0.475 -0.449                              
## DHW:SpcsS.s  0.482 -0.901 -0.379 -0.597 -0.370  0.798                       
## DHW:IntlDmn  0.540 -0.852 -0.425 -0.440 -0.596  0.754    0.701              
## SpcsO.fv:ID  0.627 -0.351 -0.862 -0.442 -0.691  0.483    0.256  0.412       
## DHW:SpO.:ID -0.414  0.653  0.507  0.337  0.456 -0.798   -0.537 -0.766 -0.597
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
  anova(bleachingmod3)#no significant difference in bleaching (DHW interaction with species) between species p=0.818834, but a significant interaction between DHW and initial symbiont type (p=0.00421), with nond-hosting corals bleaching more. 
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
## DHW                    22.5923 22.5923     1 118.404 375.4414 < 2.2e-16 ***
## Species                 0.0952  0.0476     2  56.230   0.7907  0.458514    
## InitialDom              0.0143  0.0143     1  68.629   0.2375  0.627558    
## DHW:Species             0.0241  0.0120     2 120.106   0.2002  0.818834    
## DHW:InitialDom          0.6590  0.6590     1 120.392  10.9518  0.001234 ** 
## Species:InitialDom      0.0025  0.0025     1  68.629   0.0423  0.837599    
## DHW:Species:InitialDom  0.0810  0.0810     1 120.392   1.3468  0.248139    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#now looking at any batch differences within each coral species:
  mcavbleachingmod<-glmer(Shprop~DHW*Batch+(1|Colony), data=filter(allbleach,allbleach$Species=='M.cavernosa'))
  plot_resqq(mcavbleachingmod)

  mcavbleachingmod2<-as_lmerModLmerTest(mcavbleachingmod)
  summary(mcavbleachingmod2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ DHW * Batch + (1 | Colony)
##    Data: filter(allbleach, allbleach$Species == "M.cavernosa")
## 
## REML criterion at convergence: -5.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4227 -0.0828  0.0082  0.0650  5.2633 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.00000  0.0000  
##  Residual             0.03716  0.1928  
## Number of obs: 55, groups:  Colony, 10
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  1.01595    0.05271 51.00000  19.274  < 2e-16 ***
## DHW         -0.04102    0.01307 51.00000  -3.138  0.00283 ** 
## Batch2      -0.02849    0.07227 51.00000  -0.394  0.69507    
## DHW:Batch2  -0.08481    0.01590 51.00000  -5.335 2.21e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) DHW    Batch2
## DHW        -0.682              
## Batch2     -0.729  0.497       
## DHW:Batch2  0.561 -0.822 -0.682
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
  anova(mcavbleachingmod2) #significant interaction between DHW and batch on proportion of symbionts retained, p=2.21e-06, with the october batch losing more.
## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF  F value    Pr(>F)    
## DHW       4.0937  4.0937     1    51 110.1738 2.410e-14 ***
## Batch     0.0058  0.0058     1    51   0.1554    0.6951    
## DHW:Batch 1.0577  1.0577     1    51  28.4652 2.207e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
recovmod<-glm(Shprop~Recov.Days*Species*InitialDom, data=allrecov, family='quasipoisson')
plot(recovmod)

recoveryfit= with(summary(recovmod), 1 - deviance/null.deviance)
recoveryfit #0.5469, R^2 for plotted quasipoisson model (without batch)
## [1] 0.5469235
recovmod2<-glmer(Shprop~Recov.Days*Species*InitialDom+(1|Colony), data=allrecov)
plot_resqq(recovmod2)

recovmod3<-as_lmerModLmerTest(recovmod2)
summary(recovmod3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ Recov.Days * Species * InitialDom + (1 | Colony)
##    Data: allrecov
## 
## REML criterion at convergence: 1250.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7901 -0.2597 -0.0424  0.0247 10.8057 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept)  0.1652  0.4064  
##  Residual             17.4242  4.1742  
## Number of obs: 217, groups:  Colony, 29
## 
## Fixed effects:
##                                               Estimate Std. Error        df
## (Intercept)                                    3.67046    1.79581 120.33950
## Recov.Days                                    -0.12507    0.03986 188.75632
## SpeciesO.faveolata                            -3.00438    2.22969 110.98030
## SpeciesS.siderea                              -3.44919    1.48290 124.40494
## InitialDomnond                                -3.65091    1.65218 130.04791
## Recov.Days:SpeciesO.faveolata                  0.14223    0.05070 188.65474
## Recov.Days:SpeciesS.siderea                    0.15282    0.03292 187.61232
## Recov.Days:InitialDomnond                      0.16279    0.03721 188.75877
## SpeciesO.faveolata:InitialDomnond              2.96527    2.31684 114.40581
## Recov.Days:SpeciesO.faveolata:InitialDomnond  -0.17354    0.05372 188.60438
##                                              t value Pr(>|t|)    
## (Intercept)                                    2.044  0.04315 *  
## Recov.Days                                    -3.138  0.00197 ** 
## SpeciesO.faveolata                            -1.347  0.18058    
## SpeciesS.siderea                              -2.326  0.02164 *  
## InitialDomnond                                -2.210  0.02887 *  
## Recov.Days:SpeciesO.faveolata                  2.805  0.00556 ** 
## Recov.Days:SpeciesS.siderea                    4.642 6.48e-06 ***
## Recov.Days:InitialDomnond                      4.376 2.00e-05 ***
## SpeciesO.faveolata:InitialDomnond              1.280  0.20318    
## Recov.Days:SpeciesO.faveolata:InitialDomnond  -3.231  0.00146 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rcv.Dy SpcsO. SpcsS. IntlDm Rc.D:SO. R.D:SS R.D:ID SO.:ID
## Recov.Days  -0.748                                                          
## SpecsO.fvlt -0.805  0.603                                                   
## SpecisS.sdr -0.825  0.625  0.664                                            
## InitilDmnnd -0.920  0.702  0.741  0.694                                     
## Rcv.Dys:SO.  0.588 -0.786 -0.746 -0.491 -0.552                              
## Rcv.Dys:SS.  0.624 -0.826 -0.503 -0.756 -0.543  0.649                       
## Rcv.Dys:InD  0.692 -0.933 -0.557 -0.536 -0.752  0.734    0.718              
## SpcsO.fv:ID  0.656 -0.500 -0.867 -0.495 -0.713  0.655    0.387  0.536       
## Rc.D:SO.:ID -0.479  0.646  0.643  0.371  0.521 -0.869   -0.497 -0.693 -0.746
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
anova(recovmod3)# symbiont recovery was significantly differnet between species (recovery days * symbiont density interaction), with ofav p=0.00556 and ssid p=6.48e-06 recovering more than mcav for a given amount of recovery (likely indicative of the delay in symbiont recovery seen in mcav). There was also a significant interaction between initial symbiont genus and recovery days p=2.00e-05, with corals initially not hosting Durusdinium recovering more with a given amount of recovery. 
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Recov.Days                    336.61  336.61     1 188.77 19.3184 1.844e-05 ***
## Species                        81.15   40.57     2 103.20  2.3286 0.1025280    
## InitialDom                     61.05   61.05     1 114.41  3.5035 0.0637941 .  
## Recov.Days:Species            304.78  152.39     2 188.33  8.7458 0.0002333 ***
## Recov.Days:InitialDom         139.59  139.59     1 188.60  8.0111 0.0051531 ** 
## Species:InitialDom             28.54   28.54     1 114.41  1.6381 0.2031778    
## Recov.Days:Species:InitialDom 181.85  181.85     1 188.60 10.4364 0.0014578 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    mcavrecovmod<-glmer(Shprop~Recov.Days*Batch+(1|Colony), data=filter(allrecov,allrecov$Species=='M.cavernosa'))
  plot_resqq(mcavrecovmod)

  mcavrecovmod2<-as_lmerModLmerTest(mcavrecovmod)
  summary(mcavrecovmod2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ Recov.Days * Batch + (1 | Colony)
##    Data: filter(allrecov, allrecov$Species == "M.cavernosa")
## 
## REML criterion at convergence: 416.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5973 -0.4961 -0.1152  0.1772  5.0981 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.1195   0.3457  
##  Residual             8.3119   2.8830  
## Number of obs: 82, groups:  Colony, 10
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       -0.44603    0.75590 66.15752  -0.590  0.55716    
## Recov.Days         0.08342    0.01772 70.99172   4.706 1.21e-05 ***
## Batch2             0.23704    0.98147 74.25824   0.242  0.80982    
## Recov.Days:Batch2 -0.06258    0.02139 70.90131  -2.926  0.00461 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rcv.Dy Batch2
## Recov.Days  -0.769              
## Batch2      -0.753  0.592       
## Rcv.Dys:Bt2  0.637 -0.829 -0.755
  anova(mcavrecovmod2) #there was a significant interaction p=0.004611 of batch on the relationship between symbiont density and recovery days, with the October batch recovering less.
## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Recov.Days       197.495 197.495     1 70.951 23.7606 6.437e-06 ***
## Batch              0.485   0.485     1 74.258  0.0583  0.809820    
## Recov.Days:Batch  71.158  71.158     1 70.901  8.5610  0.004611 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  Spec.labs=c('M. cavernosa','O. faveolata','S. siderea')
names(Spec.labs)=c('M.cavernosa','O.faveolata','S.siderea')

SHUFFLING

Figure 2: Inter- and intra-species variation in the shuffling towards durusdinium after recovery from heat stress.

Comparison of Proportion D before start of heat stress compared to two months into recovery. Data are binned at intervals of 0.02 (for variable Proportion D), and size aesthetic relates to number of cores in each bin to reduce overplotting. Data from April and October are both included here.

mcavshift<-filter(mcav_wide_batches,!is.na(postDcat),!is.na(preDcat))
ofavshift<-filter(ofav_wide_batches,!is.na(postDcat),!is.na(preDcat))
ssidshift<-filter(ssid_wide_batches,!is.na(postDcat),!is.na(preDcat))
beforeaftershift<-rbind(mcavshift,ofavshift,ssidshift)
beforeaftershift$Treatment<-factor(beforeaftershift$Treatment,levels=c('Manipulated','Control'))
gaindd=expression(paste("Gained",italic(" Durusdinium")))
lostd=expression(paste("Lost",italic(" Durusdinium")))

   ggplot()+geom_count(data=(beforeaftershift),
     aes(x=preDcat,y=postDcat,colour=Treatment))+
facet_grid(cols=vars(Species),labeller=labeller(Species=Spec.labs))+
geom_abline(slope=1,intercept = 0,linetype='dotted')+
  scale_color_manual(values=c('brown2','blue3'),labels=c('Bleached','Control'))+
  theme_minimal(base_size = 15)%+replace%
    theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank())+
  labs(x='Proportion *Durusdinium* before',y='Proportion *Durusdinium* after')+
 guides(size=guide_legend(title='Binned n'))+
  scale_size(range=c(1,14),breaks=(c(1,2,6,10)))+
  theme(axis.title.x = element_markdown(),axis.title.y=element_markdown())+
  theme(strip.text.x = element_text(face = 'italic'),
      panel.spacing = unit(1.5, "lines"), legend.position = 'bottom', legend.title=element_text(size=12))+
      guides(colour = guide_legend(title=''))

#ggsave('redblue_shift.pdf',device='pdf', width=10,height=5) #add 'gained durusdinium'/'lost durusdinium' labels

Figure 3a: Shuffling towards high proportions of Durusdinium after heat stress recovery was greatest in M. cavernosa, followed by O. faveolata, followed by S. siderea.

Emulating the ‘symbiont shuffling’ plot to compare coral species, as in Cunning et al 2018 figure 3b. In order to be able to include mcav, which has no variation in initial proportion d, the following mcav models are independent of initial proportion d, then integrated in with the previous ofav & ssid predicted effects.

shufflemod=glm(post_propD~pre_propD + Treatment*Species*Batch, data=filter(batches,batches$Species!='M.cavernosa'), family='quasibinomial')
#NB a mixed effects model cannot be used here, since a quasibinomial distribution is needed to bound this metric between -1 and 1. 

# have a look at the data to see data distribution for the change in proportion Durusdinium
ggplot(data=filter(batches,batches$Species!='M.cavernosa'),aes(y=post_propD, x=pre_propD))+
  geom_smooth(aes(colour=Species, linetype=Treatment), method='glm',method.args = list(family = "quasibinomial"), alpha=0.1)+
  geom_point(aes(colour=Species, shape=Treatment))+
  geom_abline(slope=1,intercept = 0,linetype='dotted')+
  theme_minimal()

summary(shufflemod)# Model summary
## 
## Call:
## glm(formula = post_propD ~ pre_propD + Treatment * Species * 
##     Batch, family = "quasibinomial", data = filter(batches, batches$Species != 
##     "M.cavernosa"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.99184  -0.11958   0.03595   0.24114   1.42314  
## 
## Coefficients:
##                                                Estimate Std. Error t value
## (Intercept)                                      0.9815     0.4578   2.144
## pre_propD                                        8.3438     1.9480   4.283
## TreatmentControl                                -4.8318     2.1137  -2.286
## SpeciesS.siderea                                -1.5428     0.6744  -2.288
## BatchOctober                                     0.1981     0.9471   0.209
## TreatmentControl:SpeciesS.siderea               -0.5660     2.1306  -0.266
## TreatmentControl:BatchOctober                   -0.8514     3.1201  -0.273
## SpeciesS.siderea:BatchOctober                   -0.6365     1.4792  -0.430
## TreatmentControl:SpeciesS.siderea:BatchOctober   3.7699     3.7441   1.007
##                                                Pr(>|t|)    
## (Intercept)                                      0.0362 *  
## pre_propD                                      7.01e-05 ***
## TreatmentControl                                 0.0259 *  
## SpeciesS.siderea                                 0.0258 *  
## BatchOctober                                     0.8350    
## TreatmentControl:SpeciesS.siderea                0.7915    
## TreatmentControl:BatchOctober                    0.7859    
## SpeciesS.siderea:BatchOctober                    0.6686    
## TreatmentControl:SpeciesS.siderea:BatchOctober   0.3182    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 0.4056584)
## 
##     Null deviance: 70.630  on 66  degrees of freedom
## Residual deviance: 24.508  on 58  degrees of freedom
##   (8 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 8
#treatment had significant effect on shuffling (p=0.0259). Batch did not significantly affect shuffling (p=0.8350).
                              
#Get fitted values averaging across initial proportion durusdinium
eff <- effect(c('pre_propD', 'Species','Treatment'), shufflemod, 
              xlevels=list(pre_propD=seq(0, 1, by=0.01)))
# Get all fitted values and subsets for each treatment
res <- droplevels(data.frame(eff))

res$Species <- factor(res$Species, levels=c("O.faveolata", "S.siderea"))
res.Bl <- subset(data.frame(eff), Treatment=="Manipulated")
res.Ct <- subset(data.frame(eff), Treatment=="Control")
# Get AUC for fitted values, lower and upper confidence limits
auc <- aggregate(res[, c("fit", "lower", "upper")], 
                 by=list(Species=res$Species, Treatment=res$Treatment), 
                 FUN=function(x) (mean(x)-0.5)/0.5)
                
auc.list <- split(auc, list(auc$Treatment))
auc$Treatment=factor(auc$Treatment, levels=c('Manipulated','Control'))
levels(auc$Treatment)<-c('Bleached','Control')
auc$Species=factor(auc$Species, levels=c('O.faveolata', 'S.siderea'))
#using model independent of initial prop d for mcav
shufflemod2=lm(post_propD ~Treatment*Batch,
                 data=subset(batches, batches$Species=='M.cavernosa'))
summary(shufflemod2)
## 
## Call:
## lm(formula = post_propD ~ Treatment * Batch, data = subset(batches, 
##     batches$Species == "M.cavernosa"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63488  0.00000  0.00129  0.08423  0.12081 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    0.87919    0.04448  19.764  < 2e-16 ***
## TreatmentControl              -0.87919    0.08671 -10.139 8.19e-12 ***
## BatchOctober                   0.11952    0.06291   1.900    0.066 .  
## TreatmentControl:BatchOctober -0.10086    0.12263  -0.822    0.417    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1664 on 34 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.8729, Adjusted R-squared:  0.8617 
## F-statistic: 77.83 on 3 and 34 DF,  p-value: 2.631e-15
#again for the initial durusdinium independent model, batch is not a significant driver of shuffling, p= 0.066.  This model did not fit a quasibinomial distrubution, but is amost entirely bond between 0 and 1 regardless due to very large effect size & minimal variation in the data. 

ggplot(data=filter(batches,batches$Species=='M.cavernosa'),aes(y=post_propD, x=Treatment))+
  geom_point(aes(shape=Treatment))+
  theme_minimal()

eff2 <- effect(c('Treatment'), shufflemod2)

# Get all fitted values and subsets for each treatment level
res2 <- droplevels(data.frame(eff2))

res2.Bl <- subset(data.frame(eff2), Treatment=="Manipulated")
res2.Ct <- subset(data.frame(eff2), Treatment=="Control")
# Get AUC for fitted values, lower and upper confidence limits
auc2 <- aggregate(res2[, c("fit", "lower", "upper")], 
                 by=list(Treatment=res2$Treatment),
                 FUN=function(x) (mean(x)))
levels(auc2$Treatment)<-c('Control','Bleached')
auc2.list <- split(auc2, list(auc2$Treatment))
auc2$Species=factor(rep('M.cavernosa'))

speciesshuff<-rbind(auc,auc2)
speciesshuff$Species=factor(speciesshuff$Species, levels=c('M.cavernosa', 'O.faveolata', 'S.siderea'))

#quote the shuffling metric for bleached cores of each species, 'S': mcav=0.938948279, ofav=0.927486547, ssid=0.722313707. 

ggplot(data=speciesshuff)+
  geom_hline(yintercept=0,linetype='dashed')+
  geom_errorbar(aes(ymin=lower, ymax=upper, x=Species, colour=Species, group=Treatment),size=0.5, position=position_dodge(width=0.5), width=0.2)+
   geom_point(aes(y=fit, x=Species,shape=Treatment, colour=Species),size=2, position=position_dodge(width=0.5))+
  theme_minimal(base_size = 13)%+replace%
   theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.minor.x = element_blank(),
         panel.grid.major.x = element_blank(),
          panel.grid.minor.y = element_blank())+
  scale_color_manual(values=c('deeppink2','darkorange1', 'darkturquoise'))+
  guides(colour=F, shape=guide_legend(title=''))+
  theme(legend.position=c(0.1,0.15))+
  labs(y='Symbiont Shuffling', x='')+
  scale_y_continuous(limits=c(-1,1.01),expand=c(0,0))+
   scale_x_discrete(labels=c(mcav,ofav,ssid), expand=expansion(mult=c(0.5,0.2)))

#ggsave('speciesshuffle.pdf',device='pdf', width=7,height=5)
#add 'more durusdinium/less durusdinium' labels post-save

shuffling and photochemical efficiency

Differnet photochemical disadvantages of hosting Durusdinium at ambient temperatures cannot explain the difference in shuffling between O. faveolata and S. siderea.

Recreating Cunning et al 2018 test of the relative photochemical disadvantages of hosting durusdinium at ambient temperatures to explain species hierarchy (ofav>ssid) in shuffling.

ipam=filter(allbleach,allbleach$Timepoint=='Pre-bleach',allbleach$Species!='M.cavernosa')

ipaminitmod=lmer(Y2~PropD*Species*Batch+(1|Colony),data=ipam)
plot_resqq(ipaminitmod)

summary(ipaminitmod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y2 ~ PropD * Species * Batch + (1 | Colony)
##    Data: ipam
## 
## REML criterion at convergence: -203.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.86465 -0.49506 -0.03011  0.50441  1.80649 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Colony   (Intercept) 0.0001145 0.01070 
##  Residual             0.0002989 0.01729 
## Number of obs: 51, groups:  Colony, 19
## 
## Fixed effects:
##                                 Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                    0.5293324  0.0070538 25.0437420  75.042  < 2e-16
## PropD                         -0.0032823  0.0134075 27.9645585  -0.245 0.808391
## SpeciesS.siderea               0.0001896  0.0108844 26.2466965   0.017 0.986235
## Batch2                         0.0363517  0.0096538 37.8622808   3.766 0.000564
## PropD:SpeciesS.siderea        -0.0216401  0.0186729 31.1020055  -1.159 0.255317
## PropD:Batch2                  -0.0017877  0.0163661 33.5505728  -0.109 0.913669
## SpeciesS.siderea:Batch2       -0.0099099  0.0158472 33.9985376  -0.625 0.535923
## PropD:SpeciesS.siderea:Batch2 -0.0053723  0.0239426 32.1793643  -0.224 0.823879
##                                  
## (Intercept)                   ***
## PropD                            
## SpeciesS.siderea                 
## Batch2                        ***
## PropD:SpeciesS.siderea           
## PropD:Batch2                     
## SpeciesS.siderea:Batch2          
## PropD:SpeciesS.siderea:Batch2    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PropD  SpcsS. Batch2 PrD:SS. PrD:B2 SS.:B2
## PropD       -0.554                                           
## SpecisS.sdr -0.648  0.359                                    
## Batch2      -0.517  0.289  0.335                             
## PrpD:SpcsS.  0.398 -0.718 -0.621 -0.207                      
## PropD:Btch2  0.317 -0.564 -0.205 -0.634  0.405               
## SpcsS.sd:B2  0.315 -0.176 -0.420 -0.609  0.246   0.386       
## PrpD:SS.:B2 -0.217  0.386  0.312  0.434 -0.526  -0.684 -0.731
#UPDATE FOR ROSS: No significant differnece in the interaction between proportion durusdinium and coral species on Fv/Fm (p=0.255).

Figure S2: Despite greater symbiont loss in the October M. cavernosa corals, shuffling between April and October were similar for all three coral species.

#before getting fitted values, test statistical significance of batch for each species inidivually:
ofavbatchmod=glm(post_propD ~ pre_propD+Batch,
                 data=filter(batches, batches$Species=='O.faveolata', Treatment=='Manipulated'), family='quasibinomial')
ssidbatchmod=glm(post_propD ~ pre_propD + Batch,
                 data=filter(batches, batches$Species=='S.siderea', Treatment=='Manipulated'),family='quasibinomial')  
mcavbatchmod=lm(post_propD ~Batch,
                 data=filter(batches, batches$Species=='M.cavernosa', Treatment=='Manipulated'))
summary(ofavbatchmod) #p=0.683
## 
## Call:
## glm(formula = post_propD ~ pre_propD + Batch, family = "quasibinomial", 
##     data = filter(batches, batches$Species == "O.faveolata", 
##         Treatment == "Manipulated"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2216   0.0000   0.0000   0.1365   1.1338  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)    0.1034     0.5202   0.199    0.845
## pre_propD    695.2104   487.7503   1.425    0.170
## BatchOctober   0.4101     0.9890   0.415    0.683
## 
## (Dispersion parameter for quasibinomial family taken to be 0.3314469)
## 
##     Null deviance: 14.570  on 21  degrees of freedom
## Residual deviance:  7.007  on 19  degrees of freedom
##   (3 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 15
summary(ssidbatchmod) #p=0.9768
## 
## Call:
## glm(formula = post_propD ~ pre_propD + Batch, family = "quasibinomial", 
##     data = filter(batches, batches$Species == "S.siderea", Treatment == 
##         "Manipulated"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6993   0.1044   0.1079   0.1944   1.3517  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -0.40083    0.53355  -0.751   0.4605  
## pre_propD     5.57625    2.20431   2.530   0.0191 *
## BatchOctober  0.03385    1.14948   0.029   0.9768  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 0.4851499)
## 
##     Null deviance: 23.231  on 24  degrees of freedom
## Residual deviance: 12.690  on 22  degrees of freedom
##   (5 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 7
summary(mcavbatchmod) #p=0.107
## 
## Call:
## lm(formula = post_propD ~ Batch, data = filter(batches, batches$Species == 
##     "M.cavernosa", Treatment == "Manipulated"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63488  0.00062  0.00129  0.11015  0.12081 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.87919    0.05068  17.348 8.19e-16 ***
## BatchOctober  0.11952    0.07167   1.668    0.107    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1896 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.09662,    Adjusted R-squared:  0.06187 
## F-statistic: 2.781 on 1 and 26 DF,  p-value: 0.1074
eff <- effect(c('pre_propD', 'Species','Treatment', 'Batch'), shufflemod, 
              xlevels=list(pre_propD=seq(0, 1, by=0.01)))
# Get all fitted values and subsets for each treatment and now also batch
res <- droplevels(data.frame(eff))

res$Species <- factor(res$Species, levels=c("O.faveolata", "S.siderea"))
res.Bl <- subset(data.frame(eff), Treatment=="Manipulated")
res.Ct <- subset(data.frame(eff), Treatment=="Control")
res.Ap<- subset(data.frame(eff), Batch=='April')
res.Oc<- subset(data.frame(eff), Batch=='October')
# Get AUC for fitted values, lower and upper confidence limits
auc <- aggregate(res[, c("fit", "lower", "upper")], 
                 by=list(Species=res$Species, Treatment=res$Treatment, Batch=res$Batch), 
                 FUN=function(x) (mean(x)-0.5)/0.5)

auc$Treatment=factor(auc$Treatment, levels=c('Manipulated','Control'))
levels(auc$Treatment)<-c('Bleached','Control')
auc$Species=factor(auc$Species, levels=c('O.faveolata', 'S.siderea'))
auc$Batch=factor(auc$Batch, levels=c('April','October'))



eff2 <- effect(c('Treatment','Batch'), shufflemod2)
## NOTE: TreatmentBatch is not a high-order term in the model
res2 <- droplevels(data.frame(eff2))
res2.Bl <- subset(data.frame(eff2), Treatment=="Manipulated")
res2.Ct <- subset(data.frame(eff2), Treatment=="Control")
res2.Ap<- subset(data.frame(eff2), Batch=='April')
res2.Oc<- subset(data.frame(eff2), Batch=='October')

# Get AUC for fitted values, lower and upper confidence limits
auc2 <- aggregate(res2[, c("fit", "lower", "upper")], 
                 by=list(Treatment=res2$Treatment, Batch=res2$Batch),
                 FUN=function(x) (mean(x)))
levels(auc2$Treatment)<-c('Control','Bleached')
auc2$Species=factor(rep('M.cavernosa'))
batchshuff<- rbind(auc2, auc)
batchshuff$spbatch<- interaction (batchshuff$Batch,batchshuff$Species)

ggplot(data=filter(batchshuff,batchshuff$Treatment=='Bleached'))+
  geom_errorbar(aes(ymin=lower, ymax=upper, x=Species, colour=Species, group=Batch),size=0.5, position=position_dodge(width=0.5), width=0.2)+
   geom_point(aes(y=fit, x=Species, shape=Batch, colour=Species),size=2, position=position_dodge(width=0.5))+
  theme_minimal(base_size = 13)%+replace%
   theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.minor.x = element_blank(),
         panel.grid.major.x = element_blank(),
          panel.grid.minor.y = element_blank())+
  scale_color_manual(values=c('deeppink2','darkorange1', 'darkturquoise'))+scale_x_discrete(labels=c(mcav,ofav,ssid))+
  guides(colour=F, shape=guide_legend(title=''))+
  theme(legend.position=c(0.1,0.15))+
  labs(y='Symbiont Shuffling in Bleached Corals', x='')+
coord_cartesian(clip='on', ylim=c(0,1.0))+
  scale_y_continuous(expand=c(0,0))

#ggsave('shufflebatches.pdf',device='pdf',width=7,height=5)

Figure 3b: Shuffling occurred later in recovery in M. cavernosa compared to O. faveolata and S. siderea.

shuffletimes=read.csv('shuffle_timepoints.csv',header =T)
shuffletimes$Treatment=factor(shuffletimes$Treatment)
shuff=read.csv('shuff.csv',header = T)
shuff$Colony=factor(shuff$Colony)
shuff$Species=factor(shuff$Species)
shuff$Timepoint=as.integer(shuff$Timepoint,length=3)

# have a look at the data to see data distribution for the change in proportion Durusdinium
ggplot(data=filter(shuff,shuff$Species!='M.cavernosa'),aes(y=PropD, x=PropD0))+
  geom_smooth(aes(colour=Species, linetype=factor(Timepoint)), method='glm',method.args = list(family = "quasibinomial"), alpha=0.1)+
  geom_point(aes(colour=Species, shape=factor(Timepoint)))+
  geom_abline(slope=1,intercept = 0,linetype='dotted')+
  theme_minimal()
## `geom_smooth()` using formula 'y ~ x'

shuffmod=glm(PropD~PropD0+Species*Timepoint,data=filter(shuff, shuff$Species!='M.cavernosa'), family='quasibinomial')
summary(shuffmod)
## 
## Call:
## glm(formula = PropD ~ PropD0 + Species * Timepoint, family = "quasibinomial", 
##     data = filter(shuff, shuff$Species != "M.cavernosa"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.93613  -0.07129   0.06986   0.56984   1.43474  
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 -1.6604     0.6125  -2.711 0.007623 ** 
## PropD0                       6.6233     1.4531   4.558 1.18e-05 ***
## SpeciesS.siderea             0.6491     1.0954   0.593 0.554508    
## Timepoint                    1.0516     0.3087   3.407 0.000877 ***
## SpeciesS.siderea:Timepoint  -0.7745     0.4916  -1.576 0.117568    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 0.4984103)
## 
##     Null deviance: 125.769  on 133  degrees of freedom
## Residual deviance:  73.193  on 129  degrees of freedom
##   (29 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 7
anova(shuffmod,test='F')
## Analysis of Deviance Table
## 
## Model: quasibinomial, link: logit
## 
## Response: PropD
## 
## Terms added sequentially (first to last)
## 
## 
##                   Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
## NULL                                133    125.769                      
## PropD0             1   43.995       132     81.774 88.2707 2.711e-16 ***
## Species            1    1.705       131     80.069  3.4219  0.066626 .  
## Timepoint          1    5.637       130     74.432 11.3095  0.001015 ** 
## Species:Timepoint  1    1.239       129     73.193  2.4863  0.117290    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
propD1=glm(PropD ~ PropD0+Timepoint*Species,
             data=filter(shuff, shuff$Species!='M.cavernosa'), family='quasibinomial')

eff1 <- effect(c('PropD0', 'Timepoint', 'Species'), propD1, 
              xlevels=list(PropD0=seq(0, 1, by=0.01)))

res <- droplevels(data.frame(eff1))
res$Species <- factor(res$Species, levels=c("O.faveolata", "S.siderea"))
res$Timepoint <- factor(res$Timepoint,levels=c('1','2','3'))
# Get AUC for fitted values, lower and upper confidence limits
auc1 <- aggregate(res[, c("fit", "lower", "upper")], 
                 by=list(Species=res$Species,Timepoint=res$Timepoint), 
                 FUN=function(x) (mean(x)-0.5)/0.5)

#now an intial prop d independent model for mcav
propD2=lm(PropD ~ Timepoint,data=filter(shuff, shuff$Species=='M.cavernosa'))

eff2 <- effect(c('Timepoint'), propD2)

res <- droplevels(data.frame(eff2))
res$Timepoint <- factor(res$Timepoint,levels=c('1','2','3'))
# Get AUC for fitted values, lower and upper confidence limits
auc2 <- aggregate(res[, c("fit", "lower", "upper")], 
                 by=list(Timepoint=res$Timepoint), 
                 FUN=function(x) (mean(x)))
auc2$Species=factor(rep('M.cavernosa'))

aucall <- rbind(auc1, auc2)
aucall$Timepoint=factor(aucall$Timepoint, levels=c(1,2,3))
aucall$Species=factor(aucall$Species, levels=c('M.cavernosa','O.faveolata','S.siderea'))

##########
ofav=expression(paste(italic("O. faveolata")))
ssid=expression(paste(italic("S. siderea")))
mcav=expression(paste(italic("M. cavernosa")))

ggplot(aucall, aes(x = Timepoint, y = fit, group = Species)) +
  geom_errorbar(data=aucall, aes(ymin = lower, ymax = upper), 
                position = position_dodge(0.2), lwd = 0.2, width = 0.2) +
  geom_point(aes(color = Species),
             position = position_dodge(0.2), size = 2.5)+
  geom_hline(yintercept = 0, lwd = 0.1) +
  scale_x_discrete(labels=c('Post heat stress','1 month recovery','2 month recovery'),expand=expansion(mult=c(0.2,0.4)))+
  theme_minimal(base_size = 13)%+replace%
   theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.minor.x = element_blank(),
         panel.grid.major.x = element_blank(),
          panel.grid.minor.y = element_blank())+
  geom_line(linetype='dashed',alpha=0.6,aes(colour=Species),position = position_dodge(0.2))+
  scale_colour_manual(values=c('deeppink2','darkorange1','darkturquoise'),labels=c(mcav,ofav,ssid), name='')+
  coord_cartesian(clip='on', ylim=c(-.2,1.0))+
  labs(y='Cumulative Symbiont Shuffling',x='')+
  annotate(geom='text',x=3.4,y=0.1,label=gaindd,size=4)+
   annotate(geom='text',x=3.4,y=-0.1,label=lostd,size=4)+
  theme(legend.position = 'bottom')+
  scale_y_continuous(expand=c(0,0))

#ggsave('shuffletimingcumulative.pdf',device='pdf',height=5,width=7)

###UPDATE FOR ROSS: The predicted values are based of a model which controls for pre-heat stress proportion durusdinium (rather than proportion durusdinium at the previous timepoint). Making 'timepoint' an integer rather than a factor in the model also changed things. Would like to perform stats on the interaction effect on shuffling between timepoint and species, but unsure how since mcav fitted to a separate model. The '2-month recovery' points are very slightly differnet from those in the bleached Vs control shuffling plot, but I think this is indicative of shuffling at each timepoint being linked to each coral (rather than averaged within each species) for these predicted values, arther than a mistake in the code. 

Figure 3c: O. faveolata and S.siderea corals that initially hosted background Durusdinium shuffled more quickly and more fully to Durusdinium than those with no initial Durusdinium, which shuffled less than M.cavernosa that contained no initial Durusdinium.

#now look at timing for switching and shuffling
shufflegroups<- filter(shuffletimes,shuffletimes$Treatment=='Manipulated')
shufflegroups$initiald<-paste(shufflegroups$PropD1)
shufflegroups<-filter(shufflegroups,shufflegroups$initiald!='NA')
shufflegroups$initiald[shufflegroups$Species=='M.cavernosa' & shufflegroups$PropD1==0]='MC0'
shufflegroups$initiald[shufflegroups$Species!='M.cavernosa' & shufflegroups$PropD1==0]='OFSS0'
shufflegroups$initiald[shufflegroups$Species!='M.cavernosa' & shufflegroups$PropD1>0]='OFSS0TO50'
shufflegroups$initiald[shufflegroups$Species!='M.cavernosa' & shufflegroups$PropD1>0.50]='OFSS50TO100'
shufflegroups$initiald=factor(shufflegroups$initiald, levels=c('MC0','OFSS0','OFSS0TO50','OFSS50TO100'))
shufflegroups<-filter(shufflegroups,shufflegroups$initiald!='NA')

 

shuffmod4=glm(PropD4~initiald,data=shufflegroups, family='quasibinomial')
summary(shuffmod4)
## 
## Call:
## glm(formula = PropD4 ~ initiald, family = "quasibinomial", data = shufflegroups)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.04154   0.07806   0.09783   0.36168   1.32585  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           2.6943     0.5735   4.698 1.27e-05 ***
## initialdOFSS0        -3.0367     0.7139  -4.253 6.40e-05 ***
## initialdOFSS0TO50    -0.6697     0.8689  -0.771    0.443    
## initialdOFSS50TO100   3.0978     2.8069   1.104    0.274    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 0.5267079)
## 
##     Null deviance: 48.819  on 73  degrees of freedom
## Residual deviance: 28.109  on 70  degrees of freedom
##   (6 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 8
anova(shuffmod4,test='F')# prop d significantly higher after mcav switching compared to ofav/ssid switching p=6.40e-05***
## Analysis of Deviance Table
## 
## Model: quasibinomial, link: logit
## 
## Response: PropD4
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
## NULL                        73     48.819                     
## initiald  3    20.71        70     28.109 13.107 6.939e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shufflegroups= pivot_longer(data=shufflegroups, cols=starts_with('PropD'),
               names_to='Timepoint', names_prefix='PropD', values_to='PropD',
               values_drop_na=T)

MCswitchex=expression(paste( italic("M. cavernosa "),'with 0.0                             '))
OFSSswitchex=expression(paste( italic("O. faveolata & S. siderea "),'with 0.0  '))
OFSSshuffex1=expression(paste( italic("O. faveolata & S. siderea "),'with >0.0 ≤ 0.5'))
OFSSshuffex2=expression(paste( italic("O. faveolata & S. siderea "),'with >0.5'))


ggplot(shufflegroups, aes(x = Timepoint, y = PropD, group=initiald)) +
  stat_summary(aes(colour=initiald),fun.data='mean_se',
             position = position_dodge(0.2), size=0.3, show.legend =F) +
  scale_x_discrete(labels=c('Pre heat stress','Post heat stress','1 month recovery','2 month recovery'))+
  stat_summary(geom='line', aes(linetype=initiald, colour=initiald),position = position_dodge(0.2))+
  theme_minimal(base_size = 15)%+replace%
   theme(panel.border = element_rect(colour='grey20',fill=NA),
          panel.grid.minor.x = element_blank(),
         panel.grid.major.x = element_blank(),
          panel.grid.minor.y = element_blank())+
  scale_linetype_manual(values=c('solid','dashed','dashed','dashed'),labels=c(MCswitchex,OFSSswitchex,OFSSshuffex1, OFSSshuffex2), name='Initial proportion *Durusdinium* groups')+
  scale_colour_manual(values=c('blue3','blue3','purple','brown2'),labels=c(MCswitchex,OFSSswitchex,OFSSshuffex1, OFSSshuffex2),name='Initial proportion *Durusdinium* groups')+
  labs(y="Proportion *Durusdinium*", x='')+
  guides(colour=guide_legend(nrow=2, byrow=TRUE, title.position ='top',title.hjust = 0.5),linetype=guide_legend(nrow=2, byrow=TRUE, title.position ='top',title.hjust = 0.5))+
  theme(legend.position = 'bottom', legend.title=element_markdown(), axis.title.y = element_markdown())

#ggsave('shuffswitch.pdf',device='pdf',height=5,width=7)

###UPDATE FOR ROSS: Due to small sample size, ofav and ssid are grouped together here (model finds no significant affect of species). Now this plot includes all corals.